Free time & chunking modeling
  • Source Code
  • Report a Bug
  1. Notebooks
  2. Model 3: Non-linear recovery
  3. Explore model predictions
  • About
  • Development notes
    • Notes
    • 2021-05-16 Meeting Notes
    • Model with sqrt scaled strength
    • Playground 1: Optim stuff
  • Notebooks
    • Data
      • View the data structure
      • Exploratory data analysis
      • Subject-level data
    • Model 1: Original model
      • Main results
      • Parameter identifiability
      • Sensitivity to tau
      • Experiment 3
    • Model 2: Include encoding time
      • Main results
    • Exploring model predictions
    • Model 3: Non-linear recovery
      • Explore model predictions
      • Basic fits
  • Function reference
    • Aggregate Data
    • Calculate the deviance of a model
    • Get data object from a file
    • get_data function
    • Inverse Logit Transformation
    • Logit Transformation
    • Calculate the overall deviance
    • Preprocesses the data
    • Serial Recall Model

On this page

  • Non-linear resource recovery model
  • Just the resource recovery
  • Model predictions
  • Model predictions (with lambda = 0.5)
  1. Notebooks
  2. Model 3: Non-linear recovery
  3. Explore model predictions

Explore model predictions

  • Show All Code
  • Hide All Code

  • View Source
Published

May 20, 2024

Modified

May 20, 2024

Non-linear resource recovery model

In Exploring model predictions I used a shiny app to understand how the basic model with linear resource recovery works, and I found that it does not predict an interaction between chunking and free time proactive benefits. I suspect that a model with a non-linear recovery will be different, because for a fixed time interval it will recover more resource the lower the resource level is. Let’s explore this. Resources are recovered according to the following equation:

\[ R(t) = R_{\text{current}} + (R_{\text{max}} - R_{\text{current}})(1 - e^{-\text{rate} \times t}) \]

where \(R_{\text{max}}\) is the maximum resource level, \(R_{\text{current}}\) is the current resource level, \(\text{rate}\) is the recovery rate, and \(t\) is the time elapsed.

Just the resource recovery

Plot the amount of resource recovered for a fixed time interval as a function of the current resource level. You can move the slider to change the recovery rate.

Code
resources <- function(R, rate, time, r_max = 1) {
  (r_max - R) * (1 - exp(-rate * time))
}

# Starting resource level
pars <- expand.grid(R = seq(0, 0.9, by = 0.1), rate = seq(0.01, 1, 0.01), time = seq(0, 6, by = 0.1))
pars$resources <- resources(pars$R, pars$rate, pars$time)

# ggplot with rate as frame and then convert via ggplotly to an interactive plot
library(ggplot2)
library(plotly)
p1 <- ggplot(pars, aes(time, resources+R, color = R, group = R, frame = rate)) +
  geom_line() +
  theme_test(base_size = 14) +
  labs(
    title = "Resource recovery",
    x = "Elapse time for recovery",
    y = "Available resources",
    color = "Initial resource level"
  )
ggplotly(p1)
Code
p2 <- ggplot(pars, aes(time, resources, color = R, group = R, frame = rate)) +
  geom_line() +
  theme_test(base_size = 14) +
  labs(
    title = "Resource recovery",
    x = "Elapse time for recovery",
    y = "Recovered resources",
    color = "Initial resource level"
  )
ggplotly(p2)


Here’s a fixed plot for rate = 0.5:


Code
pars <- expand.grid(R = seq(0, 0.9, by = 0.1), time = seq(0, 6, by = 0.1))
pars$resources <- resources(pars$R, 0.5, pars$time)

# ggplot with rate as frame and then convert via ggplotly to an interactive plot
library(ggplot2)
library(plotly)
p1 <- ggplot(pars, aes(time, resources+R, color = R, group = R)) +
  geom_line() +
  theme_test(base_size = 14) +
  labs(
    title = "Resource recovery",
    x = "Elapse time for recovery",
    y = "Available resources",
    color = "Initial resource level"
  )
ggplotly(p1)


Yes, now we can clearly see above that within the same time interval, the lower the resource level, the more resources are recovered. Let’s see if this produces the interaction we expected.

Model predictions

#| standalone: true
#| viewerHeight: 700
library(shiny)
library(shinylive)
library(ggplot2)
library(dplyr)
library(tidyr)

# define functions
pred_prob <- function(
    setsize, ISI = rep(0.5, setsize), item_in_ltm = rep(TRUE, setsize),
    prop = 0.2, prop_ltm = 0.5, tau = 0.15, gain = 25, rate = 0.1,
    r_max = 1, lambda = 1, growth = "linear", type = "response") {
  R <- r_max
  strengths <- vector("numeric", length = setsize)
  p_recall <- vector("numeric", length = setsize)
  prop_ltm <- ifelse(item_in_ltm, prop_ltm, 1)

  for (item in 1:setsize) {
    # strength of the item and recall probability
    strengths[item] <- (prop * R)^lambda

    # amount of resources consumed by the item
    r_cost <- strengths[item]^(1 / lambda) * prop_ltm[item]
    R <- R - r_cost

    # recover resources
    R <- switch(growth,
      "linear" = min(r_max, R + rate * ISI[item]),
      "asy" = R + (r_max - R) * (1 - exp(-rate * ISI[item]))
    )
  }

  if (type == "response") {
    1 / (1 + exp(-(strengths - tau) * gain))
  } else {
    strengths
  }
}

predict_resmodel <- function(object, data, group_by, type = "response", ...) {
  if (missing(group_by)) {
    pred <- pred_prob(
      setsize = nrow(data),
      ISI = data$ISI,
      item_in_ltm = data$item_in_ltm,
      prop = object["prop"],
      prop_ltm = object["prop_ltm"],
      tau = object["tau"],
      gain = object["gain"],
      rate = object["rate"],
      type = type,
      ...
    )
    return(pred)
  }

  by <- do.call(paste, c(data[, group_by], sep = "_"))
  out <- lapply(split(data, by), function(x) {
    x$pred_tmp_col295 <- predict_resmodel(object, x, type = type, ...)
    x
  })
  out <- do.call(rbind, out)
  out <- suppressMessages(dplyr::left_join(data, out))
  out$pred_tmp_col295
}

data <- expand.grid(
  chunk = c("known", "random"),
  gap = seq(500, 6000, by = 225),
  itemtype = c("SP1-3", "SP4-6", "SP7-9")
)

data$ISI <- ifelse(data$itemtype == "SP1-3", data$gap / 1000, 0.5)
data$item_in_ltm <- ifelse(data$itemtype == "SP1-3", data$chunk == "known", FALSE)
shinyApp(
  ui = fluidPage(
    titlePanel("Interactive Plot"),
    sidebarLayout(
      sidebarPanel(
        sliderInput("prop", "prop:", min = 0, max = 1, value = 0.15),
        sliderInput("prop_ltm", "prop_ltm:", min = 0, max = 1, value = 0.5),
        sliderInput("rate", "rate:", min = 0, max = 1, value = 0.25),
        sliderInput("gain", "gain:", min = 1, max = 100, value = 33),
        sliderInput("tau", "tau:", min = 0, max = 1, value = 0.12),
      ),
      mainPanel(
        plotOutput("distPlot")
      )
    )
  ),
  server = function(input, output) {
    output$distPlot <- renderPlot({
      par <- c(prop = input$prop, prop_ltm = input$prop_ltm, rate = input$rate, gain = input$gain, tau = input$tau)
      data |>
        mutate(
          Probability = predict_resmodel(par, data = data, group_by = c("gap", "chunk"), growth = "asy", lambda = 1),
          Strength = predict_resmodel(par, data = data, group_by = c("gap", "chunk"), type = "strength", growth = "asy", lambda = 1)
        ) |>
        pivot_longer(c(Probability, Strength), names_to = "type", values_to = "value") |>
        ggplot(aes(gap, value, color = chunk, group = chunk)) +
        geom_line() +
        scale_color_discrete("1st chunk LTM?") +
        facet_grid(type ~ itemtype, scales = "free") +
        theme_classic(base_size = 14)
    })
  }
)

Model predictions (with lambda = 0.5)

#| standalone: true
#| viewerHeight: 700
library(shiny)
library(shinylive)
library(ggplot2)
library(dplyr)
library(tidyr)

# define functions
pred_prob <- function(
    setsize, ISI = rep(0.5, setsize), item_in_ltm = rep(TRUE, setsize),
    prop = 0.2, prop_ltm = 0.5, tau = 0.15, gain = 25, rate = 0.1,
    r_max = 1, lambda = 1, growth = "linear", type = "response") {
  R <- r_max
  strengths <- vector("numeric", length = setsize)
  p_recall <- vector("numeric", length = setsize)
  prop_ltm <- ifelse(item_in_ltm, prop_ltm, 1)

  for (item in 1:setsize) {
    # strength of the item and recall probability
    strengths[item] <- (prop * R)^lambda

    # amount of resources consumed by the item
    r_cost <- strengths[item]^(1 / lambda) * prop_ltm[item]
    R <- R - r_cost

    # recover resources
    R <- switch(growth,
      "linear" = min(r_max, R + rate * ISI[item]),
      "asy" = R + (r_max - R) * (1 - exp(-rate * ISI[item]))
    )
  }

  if (type == "response") {
    1 / (1 + exp(-(strengths - tau) * gain))
  } else {
    strengths
  }
}

predict_resmodel <- function(object, data, group_by, type = "response", ...) {
  if (missing(group_by)) {
    pred <- pred_prob(
      setsize = nrow(data),
      ISI = data$ISI,
      item_in_ltm = data$item_in_ltm,
      prop = object["prop"],
      prop_ltm = object["prop_ltm"],
      tau = object["tau"],
      gain = object["gain"],
      rate = object["rate"],
      type = type,
      ...
    )
    return(pred)
  }

  by <- do.call(paste, c(data[, group_by], sep = "_"))
  out <- lapply(split(data, by), function(x) {
    x$pred_tmp_col295 <- predict_resmodel(object, x, type = type, ...)
    x
  })
  out <- do.call(rbind, out)
  out <- suppressMessages(dplyr::left_join(data, out))
  out$pred_tmp_col295
}

data <- expand.grid(
  chunk = c("known", "random"),
  gap = seq(500, 6000, by = 225),
  itemtype = c("SP1-3", "SP4-6", "SP7-9")
)

data$ISI <- ifelse(data$itemtype == "SP1-3", data$gap / 1000, 0.5)
data$item_in_ltm <- ifelse(data$itemtype == "SP1-3", data$chunk == "known", FALSE)
shinyApp(
  ui = fluidPage(
    titlePanel("Interactive Plot"),
    sidebarLayout(
      sidebarPanel(
        sliderInput("prop", "prop:", min = 0, max = 1, value = 0.29),
        sliderInput("prop_ltm", "prop_ltm:", min = 0, max = 1, value = 0.62),
        sliderInput("rate", "rate:", min = 0, max = 1, value = 0.25),
        sliderInput("gain", "gain:", min = 1, max = 100, value = 11),
        sliderInput("tau", "tau:", min = 0, max = 1, value = 0.41),
      ),
      mainPanel(
        plotOutput("distPlot")
      )
    )
  ),
  server = function(input, output) {
    output$distPlot <- renderPlot({
      par <- c(prop = input$prop, prop_ltm = input$prop_ltm, rate = input$rate, gain = input$gain, tau = input$tau)
      data |>
        mutate(
          Probability = predict_resmodel(par, data = data, group_by = c("gap", "chunk"), growth = "asy", lambda = 0.5),
          Strength = predict_resmodel(par, data = data, group_by = c("gap", "chunk"), type = "strength", growth = "asy", lambda = 0.5)
        ) |>
        pivot_longer(c(Probability, Strength), names_to = "type", values_to = "value") |>
        ggplot(aes(gap, value, color = chunk, group = chunk)) +
        geom_line() +
        scale_color_discrete("1st chunk LTM?") +
        facet_grid(type ~ itemtype, scales = "free") +
        theme_classic(base_size = 14)
    })
  }
)
Back to top
Exploring model predictions
Basic fits
Source Code
---
title: "Explore model predictions"
format: html
date: 05/20/2024
date-modified: last-modified
filters:
  - shinylive
---

## Non-linear resource recovery model

In [Exploring model predictions](model-explore.qmd) I used a shiny app to understand how the basic model with linear resource recovery works, and I found that it does not predict an interaction between chunking and free time proactive benefits. I suspect that a model with a non-linear recovery will be different, because for a fixed time interval it will recover more resource the lower the resource level is. Let's explore this. Resources are recovered according to the following equation:

$$
R(t) = R_{\text{current}} + (R_{\text{max}} - R_{\text{current}})(1 - e^{-\text{rate} \times t})
$$

where $R_{\text{max}}$ is the maximum resource level, $R_{\text{current}}$ is the current resource level, $\text{rate}$ is the recovery rate, and $t$ is the time elapsed.

## Just the resource recovery

Plot the amount of resource recovered for a fixed time interval as a function of the current resource level. You can move the slider to change the recovery rate.

```{r resource-recovery}
#| message: false
resources <- function(R, rate, time, r_max = 1) {
  (r_max - R) * (1 - exp(-rate * time))
}

# Starting resource level
pars <- expand.grid(R = seq(0, 0.9, by = 0.1), rate = seq(0.01, 1, 0.01), time = seq(0, 6, by = 0.1))
pars$resources <- resources(pars$R, pars$rate, pars$time)

# ggplot with rate as frame and then convert via ggplotly to an interactive plot
library(ggplot2)
library(plotly)
p1 <- ggplot(pars, aes(time, resources+R, color = R, group = R, frame = rate)) +
  geom_line() +
  theme_test(base_size = 14) +
  labs(
    title = "Resource recovery",
    x = "Elapse time for recovery",
    y = "Available resources",
    color = "Initial resource level"
  )
ggplotly(p1)

p2 <- ggplot(pars, aes(time, resources, color = R, group = R, frame = rate)) +
  geom_line() +
  theme_test(base_size = 14) +
  labs(
    title = "Resource recovery",
    x = "Elapse time for recovery",
    y = "Recovered resources",
    color = "Initial resource level"
  )
ggplotly(p2)
```

<br>

Here's a fixed plot for rate = 0.5:

<br>

```{r resource-recovery-fixed}
#| message: false

pars <- expand.grid(R = seq(0, 0.9, by = 0.1), time = seq(0, 6, by = 0.1))
pars$resources <- resources(pars$R, 0.5, pars$time)

# ggplot with rate as frame and then convert via ggplotly to an interactive plot
library(ggplot2)
library(plotly)
p1 <- ggplot(pars, aes(time, resources+R, color = R, group = R)) +
  geom_line() +
  theme_test(base_size = 14) +
  labs(
    title = "Resource recovery",
    x = "Elapse time for recovery",
    y = "Available resources",
    color = "Initial resource level"
  )
ggplotly(p1)
```

<br>

Yes, now we can clearly see above that within the same time interval, the lower the resource level, the more resources are recovered. Let's see if this produces the interaction we expected.

## Model predictions

::: {.column-page-right}

```{shinylive-r}
#| standalone: true
#| viewerHeight: 700
library(shiny)
library(shinylive)
library(ggplot2)
library(dplyr)
library(tidyr)

# define functions
pred_prob <- function(
    setsize, ISI = rep(0.5, setsize), item_in_ltm = rep(TRUE, setsize),
    prop = 0.2, prop_ltm = 0.5, tau = 0.15, gain = 25, rate = 0.1,
    r_max = 1, lambda = 1, growth = "linear", type = "response") {
  R <- r_max
  strengths <- vector("numeric", length = setsize)
  p_recall <- vector("numeric", length = setsize)
  prop_ltm <- ifelse(item_in_ltm, prop_ltm, 1)

  for (item in 1:setsize) {
    # strength of the item and recall probability
    strengths[item] <- (prop * R)^lambda

    # amount of resources consumed by the item
    r_cost <- strengths[item]^(1 / lambda) * prop_ltm[item]
    R <- R - r_cost

    # recover resources
    R <- switch(growth,
      "linear" = min(r_max, R + rate * ISI[item]),
      "asy" = R + (r_max - R) * (1 - exp(-rate * ISI[item]))
    )
  }

  if (type == "response") {
    1 / (1 + exp(-(strengths - tau) * gain))
  } else {
    strengths
  }
}

predict_resmodel <- function(object, data, group_by, type = "response", ...) {
  if (missing(group_by)) {
    pred <- pred_prob(
      setsize = nrow(data),
      ISI = data$ISI,
      item_in_ltm = data$item_in_ltm,
      prop = object["prop"],
      prop_ltm = object["prop_ltm"],
      tau = object["tau"],
      gain = object["gain"],
      rate = object["rate"],
      type = type,
      ...
    )
    return(pred)
  }

  by <- do.call(paste, c(data[, group_by], sep = "_"))
  out <- lapply(split(data, by), function(x) {
    x$pred_tmp_col295 <- predict_resmodel(object, x, type = type, ...)
    x
  })
  out <- do.call(rbind, out)
  out <- suppressMessages(dplyr::left_join(data, out))
  out$pred_tmp_col295
}

data <- expand.grid(
  chunk = c("known", "random"),
  gap = seq(500, 6000, by = 225),
  itemtype = c("SP1-3", "SP4-6", "SP7-9")
)

data$ISI <- ifelse(data$itemtype == "SP1-3", data$gap / 1000, 0.5)
data$item_in_ltm <- ifelse(data$itemtype == "SP1-3", data$chunk == "known", FALSE)
shinyApp(
  ui = fluidPage(
    titlePanel("Interactive Plot"),
    sidebarLayout(
      sidebarPanel(
        sliderInput("prop", "prop:", min = 0, max = 1, value = 0.15),
        sliderInput("prop_ltm", "prop_ltm:", min = 0, max = 1, value = 0.5),
        sliderInput("rate", "rate:", min = 0, max = 1, value = 0.25),
        sliderInput("gain", "gain:", min = 1, max = 100, value = 33),
        sliderInput("tau", "tau:", min = 0, max = 1, value = 0.12),
      ),
      mainPanel(
        plotOutput("distPlot")
      )
    )
  ),
  server = function(input, output) {
    output$distPlot <- renderPlot({
      par <- c(prop = input$prop, prop_ltm = input$prop_ltm, rate = input$rate, gain = input$gain, tau = input$tau)
      data |>
        mutate(
          Probability = predict_resmodel(par, data = data, group_by = c("gap", "chunk"), growth = "asy", lambda = 1),
          Strength = predict_resmodel(par, data = data, group_by = c("gap", "chunk"), type = "strength", growth = "asy", lambda = 1)
        ) |>
        pivot_longer(c(Probability, Strength), names_to = "type", values_to = "value") |>
        ggplot(aes(gap, value, color = chunk, group = chunk)) +
        geom_line() +
        scale_color_discrete("1st chunk LTM?") +
        facet_grid(type ~ itemtype, scales = "free") +
        theme_classic(base_size = 14)
    })
  }
)
```

:::


## Model predictions (with lambda = 0.5)

::: {.column-page-right}

```{shinylive-r}
#| standalone: true
#| viewerHeight: 700
library(shiny)
library(shinylive)
library(ggplot2)
library(dplyr)
library(tidyr)

# define functions
pred_prob <- function(
    setsize, ISI = rep(0.5, setsize), item_in_ltm = rep(TRUE, setsize),
    prop = 0.2, prop_ltm = 0.5, tau = 0.15, gain = 25, rate = 0.1,
    r_max = 1, lambda = 1, growth = "linear", type = "response") {
  R <- r_max
  strengths <- vector("numeric", length = setsize)
  p_recall <- vector("numeric", length = setsize)
  prop_ltm <- ifelse(item_in_ltm, prop_ltm, 1)

  for (item in 1:setsize) {
    # strength of the item and recall probability
    strengths[item] <- (prop * R)^lambda

    # amount of resources consumed by the item
    r_cost <- strengths[item]^(1 / lambda) * prop_ltm[item]
    R <- R - r_cost

    # recover resources
    R <- switch(growth,
      "linear" = min(r_max, R + rate * ISI[item]),
      "asy" = R + (r_max - R) * (1 - exp(-rate * ISI[item]))
    )
  }

  if (type == "response") {
    1 / (1 + exp(-(strengths - tau) * gain))
  } else {
    strengths
  }
}

predict_resmodel <- function(object, data, group_by, type = "response", ...) {
  if (missing(group_by)) {
    pred <- pred_prob(
      setsize = nrow(data),
      ISI = data$ISI,
      item_in_ltm = data$item_in_ltm,
      prop = object["prop"],
      prop_ltm = object["prop_ltm"],
      tau = object["tau"],
      gain = object["gain"],
      rate = object["rate"],
      type = type,
      ...
    )
    return(pred)
  }

  by <- do.call(paste, c(data[, group_by], sep = "_"))
  out <- lapply(split(data, by), function(x) {
    x$pred_tmp_col295 <- predict_resmodel(object, x, type = type, ...)
    x
  })
  out <- do.call(rbind, out)
  out <- suppressMessages(dplyr::left_join(data, out))
  out$pred_tmp_col295
}

data <- expand.grid(
  chunk = c("known", "random"),
  gap = seq(500, 6000, by = 225),
  itemtype = c("SP1-3", "SP4-6", "SP7-9")
)

data$ISI <- ifelse(data$itemtype == "SP1-3", data$gap / 1000, 0.5)
data$item_in_ltm <- ifelse(data$itemtype == "SP1-3", data$chunk == "known", FALSE)
shinyApp(
  ui = fluidPage(
    titlePanel("Interactive Plot"),
    sidebarLayout(
      sidebarPanel(
        sliderInput("prop", "prop:", min = 0, max = 1, value = 0.29),
        sliderInput("prop_ltm", "prop_ltm:", min = 0, max = 1, value = 0.62),
        sliderInput("rate", "rate:", min = 0, max = 1, value = 0.25),
        sliderInput("gain", "gain:", min = 1, max = 100, value = 11),
        sliderInput("tau", "tau:", min = 0, max = 1, value = 0.41),
      ),
      mainPanel(
        plotOutput("distPlot")
      )
    )
  ),
  server = function(input, output) {
    output$distPlot <- renderPlot({
      par <- c(prop = input$prop, prop_ltm = input$prop_ltm, rate = input$rate, gain = input$gain, tau = input$tau)
      data |>
        mutate(
          Probability = predict_resmodel(par, data = data, group_by = c("gap", "chunk"), growth = "asy", lambda = 0.5),
          Strength = predict_resmodel(par, data = data, group_by = c("gap", "chunk"), type = "strength", growth = "asy", lambda = 0.5)
        ) |>
        pivot_longer(c(Probability, Strength), names_to = "type", values_to = "value") |>
        ggplot(aes(gap, value, color = chunk, group = chunk)) +
        geom_line() +
        scale_color_discrete("1st chunk LTM?") +
        facet_grid(type ~ itemtype, scales = "free") +
        theme_classic(base_size = 14)
    })
  }
)
```

:::